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ABSTRACT 

We study the representation of solutions of the three-dimensional quasigeostrophic (QG) equations using 
Galerkin series with standard vertical modes, with particular attention to the incorporation of active surface 
buoyancy dynamics. We extend two existi ng Ga l erkin approaches (A and B) and develop a new Galerkin ap¬ 
proximation (C). Approximation A, due to lFlierlHl978h . represents the streamfunction as a truncated Galerkin 
series and defines the potential vo rticity (PV) that satisfies the inversion problem exactly. Approximation B, 
due to lTulloch and Smitfl j2009dl . represents the PV as a truncated Galerkin series and calculates the stream- 
function that satisfies the inversion problem exactly. Approximation C, the true Galerkin approximation for 
the QG equations, represents both streamfunction and PV as truncated Galerkin series, but does not satisfy 
the inversion equation exactly. The three approximations are fundamentally different unless the boundaries 
are isopycnal surfaces. We discuss the advantages and limitations of approximations A, B, and C in terms 
of mathematical rigor and conservation laws, and illustrate their relative efficiency by solving linear stability 
problems with nonzero surface buoyancy. With moderate number of modes, B and C have have superior ac¬ 
curacy than A at high wavenumbers. Because B lacks conservation of energy, we recommend approximation 
C for constructing solutions to the surface-active QG equations using Galerkin series with standard vertical 
modes. 


1. Introduction 

Recent interest in upper-ocean dynamics and sub- 
mesoscale turbulence has focussed attention on surface 
geostrophic dynamics and the role of surface buoyancy 
variations. A main issue is the representation of active sur¬ 
face buoyancy by finite vertical truncations of the quasi¬ 
equations. Standard multi-layer (e.g., 
tid modal approximations (e.g., iFlierll 
there is no variation of buoyancy on the 

surfaces. 

Here we explore the representation of surface and inte¬ 
rior dynamics using the simple and familiar vertical modes 
of physical oceanography. These modes, denoted here by 
Pn{z), are defined by the Sturm-Liouville eigenproblem 


geostrophic (QG) 


1978h assume that 


d fo dpn 
dz N^ dz 


-K;fp„, 


( 1 ) 
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with homogeneous Neumann boundary conditions at the 
bottom (z = z~) and top (z = z”^) surfaces of the domain: 

^(z±)=0. (2) 

In ([TJ A is the buoyancy frequency and /q is the Coriolis 
parameter. The eigenvalue K„ in ([T]i is the deformation 
wavenumber of the «’th mode. With normalization, the 
modes satisfy the orthogonality condition 

^ J PnPmdz — ^mn ; ( 3 ) 

where h z+ — is the depth. The barotropic mode is 
po = 1 and Kq = 0. 


The modes defined by the eigenproblem H]) and ^ pro¬ 
vide a fundamental basis for representing solutions of both 
the primitive and quas igeostroph i c equations as a linear 


combination of |d„I (Gill 

19821 Pedloskv 

19871; IVallisI 

l2006t [Ferrari and Wunschl 

ZOlOilLaCasce 2012 

). In fact. 


the set {p„} is mathematically complete and can be used 
to represent any field with finite square integral, 
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Even if the field 0 has nonzero derivative at or internal 
discontinuities, its representation as a linear combination 
of the basis functions {p„} converges inL^(z^,z+) i.e., the 
integral of the squared error goe s to zero as the number of 
basis functions increases (e.g., [Hunter and Nachtergaelel 

I 2 OOIL ch. 10). 


Despite the rigorous assurance of completeness in the 
previous paragraph, the utility of {p„} for problems with 
nonuniform surface b uoyancy has be e n questioned by sev¬ 
eral authors (e.g., L^eyr^ 120091 : [Roullet et al.l 120121 : 
ISmith and Vannestell2013l) . These authors argue that the 
homogeneous boundary conditions in @ are incompatible 
with non-zero surface buoyancy and that representation of 
the streamfunction i/r as a linear combination of {p,,} is 
useless if is non-zero on the surfaces. This supposed 
incompatibility of (|2]i with non-zero surface buoyancy is 
a main motivation for a ne w, alternative set of orthogon al 
basis functions proposed bv ISmith and Vannestel(l2013h . 

The aim of this paper is to obtain a good Galerkin ap¬ 
proximation to solutions of the QG equation with non¬ 
zero surface buoyancy using the familiar basis {p„}. We 
show that that both the inversion problem and evolution¬ 
ary dynamics can be handled using {p„} to represent the 
streamfunction. As part of this program w e revisit and 
extend two existing moda l approximations (iFlier l[T97i 
iTulloch and Smith! l2009al) . and develop a new Galerkin 
approximation. We discuss the relative merit of the 
three approximations in terms of their mathematical rigor 
and conservation laws, and illustrate their efficiency and 
caveats by solving linear stability problems with nonzero 
surface buoyancy. 

Using concrete examples, we show that the concerns 
expressed by earlier authors regarding the suitability of 
the standard modes {p„} are over-stated: even with non¬ 
zero surface buoyancy, the Galerkin expansion for xj/ in 
terms of {p„} converges absolutely and uniformly with no 
Gibbs phenomena. A modest number of terms provides a 
good approximation to i/r thi'oughout the domain, includ¬ 
ing on the top and bottom boundaries. In other words, the 
surface streamfunction can be expanded in terms of {p„} 
and, with enough modes, this representation can then be 
used to accurately calculate the advection of non-zero sur¬ 
face buoyancy. In section|5]we illustrate this procedure by 
solving the classic Eady problem using the basis {p„} for 
the streamfunction. 


2. The exact system 

In this section we summarize the basi c properties of th e 
QG system. For a detailed derivation see |Pedloskvl(ll987ll . 


a. Formulation 

The streamfunction is denoted \f/{x,y,z,t) and we use 
the following notation. 

u = -Wy, v = V4, (5) 

The variable tJ is related to the buoyancy by = N^'d-ffo. 
The QG potential vorticity (QGPV) equation is 

^tq + ^{y/,q)+Pv = 0, (6) 

where the potential vorticity is 

q={A + L)xif, (7) 


with 


+ and (8) 

Also in ([6l), the Jacobian is J(A,B) =* d^A dyB — dyA dji. 

The boundary conditions at the top (z = z+) and bottom 
(z = z^) are that w = 0, or equivalently 

= J(v/*,t7=^)=0. (9) 

Above we have used the superscripts -|- and — to denote 
evaluation atz^ andz^ e.g., 1 //+ = \\r{x,y,z^d)- 


b. Quadratic conservation laws 

In the absence of sources and sinks, the exact QG sys¬ 
tem has four quadratic conservation laws: energy, poten¬ 
tial enstrophy, and surface buoyancy varian ce at the two 
surfaces (e.g., IPedloskvIl 19871 : Ivallisll200tj) . Throughout 
we assume horizontal periodic boundary conditions. 

The well-known energy conservation law is 


where 



( 10 ) 


( 11 ) 


The total energy is poE, where po is a reference density. 
An alternative expression for E is 


E = -\j\ifqAV + \j\if^d+-\F^d-<lS. (12) 

If ^ = 0 (e.g., as in the Eady problem) then (fT2l) expresses 
E in terms of surface contributions. 

If j3 = 0 then there are many quadratic potential enstro¬ 
phy invariants: the volume integral of q^A{z), with A(z) an 
arbitrary function of the vertical coordinate, is conserved. 
The choice A(z) = 5(z — z*) reduces to conservation of the 
surface integral of q^ at any level z*. 
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IChamevI(Il97lh noted that, in a doubly periodic domain, 
nonzero j3 destroys all these quadratic potential enstrophy 
conservation laws, including the conservation of poten¬ 
tial enstrophy defined simply as the volume integral of q^. 
Multiplying the QGPV equation (|6]l by q, and integrating 
by parts, we obtain 

A I iq^dV + P I [vt7]^!d5 = 0. (13) 

The potential enstrophy equation (fTSl i is the finite-depth 
analog of Chamey’s equation (13). To make progress 
Chai'ney assumed = 0 at the ground. But the j3-term on 
the right of (fPl) can be eliminated by cross-multiplying 
the QGPV equation (|6ll evaluated at the surfaces with 
the boundary conditions (|9ll, and combining with (fT3l i. 
Thus nonzero j3 selects a uniquely conserved potential en¬ 
strophy from the infinitude of j3 = 0 potential enstrophy 
conservation laws: 


where the potential enstrophy is 

Z = j y^dV-j q+^+ -q-i}-dS. (15) 

With j3 7 ^ 0 the surface contributions in (fTSl l are required to 
form a conserved quadratic quantity involving q^. Notice 
that ( fTSl ) is not sign-definite. To our knowledge, the con¬ 
servation law in ( fldl i and (fTsI) is previously unremarked. 

Finally, in addition to E and Z, the surface buoyancy 
variance is conserved on each surface 

A/i(z7±)2d5 = 0. (16) 

Thus, with j3 7 ^ 0, the QG model has four quadratic con¬ 
servation laws: E, Z and the buoyancy variance at the two 
surfaces. 

3. Galerkin approximation using standard vertical 
modes 

A straightforward approach is to represent the stream- 
function by linearly combining the first N -f 1 vertical 
modes. The mean square error in this approximation is 


where the coefficients in the sum above are 

fn{x,y,t)=};j^ \l/Pndz. (19) 

Throughout we use the superscript" to denote a Galerkin 
coefficient defined via projection of a field onto a vertical 
mode. 

In complete analogy with the streamfunction, one can 
also develop an (N -|- l)-mode Galerkin approximation to 
the PV: 

def 

qN(.x,y,z,t) = Y,qn{x,y,t)p„{z), (20) 

n=0 

with coefficients 

fZ^ 

qn = ]lj qPndz. (21) 

The construction of the Galerkin approximation q^ above 
minimizes a mean square error err^ defined in analogy 
with (fTTl i. 

Now recall that the exact Xj/ and q are related by the 
elliptic “inversion problem”: 

iA + L)xi/ = q, (22) 

with boundary conditions at 

(#)Vz = t7±. (23) 

The Galerkin approximations in (fTSl) through (|2TI) are de¬ 
fined independently of the information in (l22l) and ( |2^ . 
The relationship between the Galerkin coefficients q„ and 
\j/„ is obtained by multiplying ( l22l i by ^Pn(z) and integrat¬ 
ing over the depth. Noting the intermediate result 

i ^ p„Lvrdz = i [p+- p- 1?-] - , (24) 

we obtain 

^„ = A„vf„-fi(p+T7+-p^i7^), (25) 

-V-' 

surface terms 


err^(flO,fli,---«N) 



(17) 


We use a roman font, and context, to distinguish the tmn- 
cation index N in (fTTI) from the buoyancy frequency N{z). 
The coefficients ao through on are determined to minimize 
sn^r, and thus one obtains the Galerkin approximation 
to the exact streamfunction: 


¥°{x,y,z,t) = £ \j/„{x,y,t)pn{z ), (18) 

n=0 


where is the n'ih mode Helmholtz operator 

A„ = A-k-2. (26) 

The relation in ( l25l l is the key to a good Galerkin approxi¬ 
mation to surface-active quasigeostrophic dyamics. 

Term-by-term differentiation of the y/^-series in (fTSl) 
does not give the q^ series in (l20l i unless = 0. In other 
words, term-by-term differentiation does not produce the 
correct relation (l25l l between q^ and \j/„. Thus the Galerkin 
tmncated PV and the Galerkin truncated streamfunction 
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do not satisfy the inversion boundary value problem ex¬ 
actly 


(A -I- L) 7^ • (27) 


Despite dZTl) . the tmncated series and are the best 
least-squares approximations to yr and q. 

Notice that, in analogy with the Galerkin approxima¬ 
tions for q and yr, 

K = and 5 - = ip-, (28) 


where 


= L P" = L p« > 

«=0 rt=0 

(29) 

are finite approximations to distributions 5(z — z^) at the 
surfaces. Of course, these surface 5-distributions do not 
satisfy the iJ convergence condition in (|4]i and thus the 
series in (|2^ only converge in a distributional sense.(e.g., 
[Hunter and Nachtergaelel2doIl) . For instance, if 0 satisfies 
the !?■ convergence condition in Q, then 


<l>{z)d^^{z)dz j 0(z)5(z-z+)dz= (/)(z+), 

(30) 

as N ^ 00 . Thus, in that limit. 


(A +L)^, (31) 


where ^ denotes distributional convergence. The right- 
hand-side ofj^lJstheBrethertonian modified potential 
vorticity (IBrethertonll 1966h with the boundary conditions 
incorporated as PV sheets. To illustrate ( l27l i and (l3Tl i we 
present an elementary example that is is relevant to our 
discussion of the Eady problem in section|5] 


An elementary example: the Eady basic state 

As an example, consider the case with constant buoy¬ 
ancy frequency N. We use nondimensional units so that 
the surfaces are at z- = — 1 and z^ = 0. The standard ver¬ 
tical modes are po = 1 and, for « > 1 

p„ = \/2cos(n7rz), (32) 


with K„ = me. 

We consider the basic state of the 
streamfunction 

V7 = -(l-fz)y, 
u 

and zero interior PV q = 0 and j3 = 0. The surface buoy¬ 
ancies are 1 }^ = —y. 

The Galerkin expansion of the PV ^ = 0 is exact: = 0 

and therefore q^ = 0. The truncated Galerkin expansion 


Eady problem with 
(33) 


of rjf follows from either ( IT9l l or ( l25l l and is 


V^N = 


jPo- 


- 2 V 2 


Pi P3 

{3n)^ 


Pn 

(N;r)" 


y- 


(34) 

(We assume that N is odd, so that the last term in the trun¬ 
cated series is as above.) Despite the nonzero derivative 
of Iff at the boundaries, the series in dMl is absolutely and 
uniformly convergent on the closed interval — 1 < z < 0. 
The N-^ behavior of the series (l34ll ensures uniform con¬ 
verge nce, e.g., using the M-test ( Hunter and Nachtergaelel 
l200lh . There are no Gibbs oscillations and a modest num¬ 
ber of terms provides a good approximation to the base 
velocity U (Eigure[T^). 

Now, to illustrate ( |27] ) and ( l3Tl i, notice that 


(A -|- L) — 2\/2(pi + P3 H- Pn)}' 

_2sin[(N-|-l)7rz] 
sin(7rz) 


-y- 


(35) 

(36) 


The series ([33 does not converge in a point- 
wise sense. Howeve r, in a distributional sense 
dHunter and Nachtergaelel I 2 OOIL ch. 11), the exact 
sum in ( l36l l does converge to 5-distributions on the 
boundaries; see figures [TJb) andfljc). These b oundary d- 
distrib utions are the Brethertonian PV sheets dBrethertonl 

[T^ . 


4. Three approximations 

In d27]) we noted that the Galerkin approximations to 
\jf and q do not exactly satisfy the inversion relation. To 
address this error there are at least three different approxi¬ 
mations one can make: 

Approximation A: Use the truncated sum in (fTSl) 
as a least-squares Galerkin approximation to the stream- 
function \f/. But do not use the Galerkin approximation 
for q. Instead, define the approximate PV, q^{x,y,z,t), so 
that the interior inversion relation is satisfied exactly: 

4"^f(A + L)v/;?. (37) 

This is the approximation introduced by Iplierll dl978h . 
which is now regarded as the standard in physical 
oceanography. Note that in (I 37 I 1 is not the least-squares 
approximation to the exact q. Moreover, the approxima¬ 
tion approaches the Brethertonian PV on the right of 
dSB as N ^ 00 . 

Approximation B: Use the truncated sum in d20l i as 
a least-squares Galerkin approximation to the PV q. But 
do not use the Galerkin approximation for y/. Instead, de¬ 
fine the approximate streamfunction, y/^{x,y,z,t), as the 
solution to the inversion boundary value problem 

(A -I- L) = 9n I 


(38) 
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Fig. 1. Nondimensional base-state for the Eady problem using various truncation for the series o. In the middle panel N is the number of 
baroclinic modes, (a) Zonal velocity: although the truncation has zero slope at the boundaries there are no Gibbs oscillations, (b) Meridional PV 
gradient associated with the truncated series (36). (c) as in (b) but with an expanded abscissa. As N increases, the PV gradient distributionally 
converges to two Brethertonian delta functions at the boundaiies. 


with boundary conditions 

= (39) 

This is the approx imation introduced by 
IXulloch and Smith! (12009b[) . Notice that (l38l l and 

(l39l l is an approximation to the exact inversion problem 
because the interior source is q^, rather than q. In other 
words, yr® is an exact solution to an approximate version 
of the inversion problem. But is not a least-squares 
approximation to the exact \jf, and nor can y/^ be written 
as a finite sum of vertical modes. 

Approximation C: Use truncated Galerkin approxima¬ 
tions and q^ for both y/ and q. In this case, as indicated 
in ( l27l i. the inversion equation will not be satisfied exactly 
by the approximate streamfunction and PV. But instead, 
one will have true least-squares approximations to both 
yf and q. To our knowledge approximation C, correctly 
accounting for the surface-buoyancy boundary terms, has 
not been previously investigated. 

In approximation A there are N + I modal amplitudes. 
In approximations B and C there are N -|- 3 degrees of free¬ 
dom: the N + I modal amplitudes q„ and the two surface 
buoyancy fields tJ*. The three approximations are equiv¬ 
alent when t?* = 0. Approximation C is a tme Galerkin 
approximation and it is important to understand its limita¬ 
tions and advantages relative to the better known alterna¬ 
tives A and B. 

Once an approximation has been chosen, one needs to 
constmct evolution equations for the Galerkin coefficients 


using the QG equations (|6]l and (|9]l. In the next three 
sub-sections, we derive evolution equations and the as¬ 
sociated inviscid conservation laws for the three approx¬ 
imations outlined above. After testing, we recommend C 
as the most reliable approximation using standard vertical 
modes. 

a. Approximation A 

Following iFlier 1 (Il978h . in approximation A the N- 
mode approximate PV is defined via ( iTTl i and, using the 
modal representation for y/§ in ( ITSl l. this is equivalent to 

. N 

^ = Z ^nfn{x,y,t) Pn{z) , (40) 

n=0 

where A„ is the Helmhol tz operator in (l26l l. Following the 
appendix of lFlierll (Il978h . one can use Galerkin projection 
of the nonlinear evolution equation (|6]l onto the modes p„ 
to obtain N + 1 evolution equations for the coefficients y/„: 

N N 

•^nmsA l^d^y/n — ^ 7 (41) 

m=0s=0 

where 



Note that cannot be computed exactly except in cases 
with simple buoyancy frequency profiles. But it suffices 
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to compute to high accuracy, e.g. using Gaussian 
q uadrature. 

Iplierll (Il978h implicitly assumed that t9^+ = tJ^ = 0, 
so that the surface terms in ( l25l l vanish and then there is 
no difference between and q^. But in general, with 
nonzero surface buoyancy, we can append evolution equa¬ 
tions for t?+ and tJ^ to approximation A. That is, in addi¬ 
tion to the N -f 1 modal equations in (IdTT i. we also have 

N 

5,t?±+Xp±J(v/„,t?±)=0. (43) 

n=0 

Above we have evaluated the i/r-series (fTSl) at z* to ap¬ 
proximate \j/^ in the surface boundary conditions. This 
approach is not satisfactory because the resulting surface 
buoyancy equations ( l43l l are dynamically passive i.e., tJ+ 
and Td-^ do not affect the interior evolution equations in 

(ED. 

Quadratic conservation laws 

Appendix A shows that Approximation A has the en¬ 
ergy conservation 

/ H^fnf+^^^fndS = 0. (44) 

n=0 

To obtain the energy analogous to E in (HD, the modal 
sum above is multiplied by the depth h. With j3 7 ^ 0, ap¬ 
proximation A has the potential enstrophy conservation 
law, 

d N „ 

-Y H^nfnYdS = 0. (45) 

M=0 

With j3 7 ^ 0, the analog of the exact potential enstrophy 
(Hsii is not conserved. Finally, with the surface equations 
in ( l43l l. approximation A also conserves surface buoyancy 
variance as in (HSll. 

b. Approximation B 

Approximation B begins with the observation that the 
exact solution of the inversion problem in (l 22 li and (| 2 D 
can be decomposed as 

\l/=<j) + o (46) 

where ^(x,y,z,t) is the “interior streamfunction” 

and g(jc,v,z,r) is t he “surface streamfunc tion” 

dLanevre and Kleiiill2006l : iTulloch and Smithll2009bh . 

The surface streamfunction o{x,y,z,t) is defined as the 
solution of the boundary value problem 

(A-fL)a = 0, (47) 

with inhomogeneous Neumann boundary conditions 

(^)"5,a(z±) = tJ±. (48) 


In approximation B, one must always solve for the surface 
streamfunction using methods other than a truncated se¬ 
ries. The solution of the surface problem (l47l i and (1481) in 
terms of standard vertical modes is 

o'= and CT,, = ^ / p„adz. (49) 

«=0 

where 

A„<7„ = -i(p+t?+-p;t?-), (50) 

and A„ is the «’th mode Helmholtz operator defined in 
(l26l) . In (HD and (l50l) we have a solution for the surface 
streamfunction a, with nonzero vertical derivative at 
the surfaces, in terms of vertical modes with zero deriva¬ 
tive. Truncations of the series ddD behave similarly to the 
tmncated series ( l34l i: convergence is absolute and uniform 
(see appendix A). 

The interior streamfunction ^{x,y,z,t) is defined as the 
solution of the boundary value problem 

(A + L)(/)=^, (51) 

with homogeneous Neumann boundary conditions 

(#)"5,0(z±)=O. (52) 

Approximation B assumes that one can solve the sur¬ 
face problem in (l47l i and (l48l l without resorting to trun¬ 
cated versions of the series in (ED. For instance, with 
constant or exponential stratific ations one can find closed- 
form, exact exp ressions for o dTulloch and Smithll2009bl : 
iLaCascel 120121) . In particular, approximation B requires 
that the two unknown Dirichlet boundary-condition func¬ 
tions a* = cr(z^) can be obtained efficiently from speci¬ 
fied Neumann boundary-condition functions and . 

The Eady problem, discussed below in section |5] is a 
prime example in which one can obtain this Neumann-to- 
Dirichlet map. 

Once a is in hand, the approximate streamfunction is 

Vn = '/’n + ^ j (53) 

where (j>§{x,y,z,t) is obtained by solving the interior in¬ 
version problem dsD with the right hand side replaced by 
the Galerkin approximation defined in (l20l i and (|2D . 
The exact solution of this approximation to the interior in¬ 
version problem is 

N 

0N = 12(54) 

n=0 

where 

d f 

= t; Pn<j>dz, and A„(^„=§„ 


(55) 
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To obtain the approximation B evolution equations we 
introduce the streamfunction (l53l l into the QGPV equation 
(|6ll and project onto mode n to obtain 

N N 

df^n^n ‘^nms J “h “h ^n) 

ni=Qs=0 

+ 121 /_P«P.vJ(ct,Aj4)dz = 0, (56) 

i =0 

with defined in (l42li . Approximation B assumes that 
the remaining integral on the second line of (l5^ can be 
evaluated exactly. This is only possible for particular mod¬ 
els of the N{z) (e.g., constant buoyancy-frequency pro¬ 
files). In practice, however, it may suffice to compute the 
integral on the second line (l56l l very accurately, e.g. using 
Gaussian quadrature. 

The evolution equations for approximation B are com¬ 
pleted with the addition of buoyancy-advection at the sur¬ 
faces 

a,tJ± + J(a±,t?±)+£ p±J(0„,tJ±)=O. (57) 

«=0 

With (l5^ and ( l57l i we have N -|- 3 evolution equations for 
the N-h3 fields 

Quadratic conservation laws 

Approximation B conserves surface buoyancy variance. 
But the conservation laws for energy and potential enstro- 
phy are problematic. The analog of the exact total energy 
(fTTl i is not generally conserved in approximation B (Ap¬ 
pendix A). With P = 0, approximation B has a potential 
enstrophy conservation 

d’^T 

-Y n{Aj„fdS = 0. (58) 

^ n=0-' 

But with j3 7 ^ 0 the analog of the exact potential enstrophy 
O is not conserved (Appendix A). 

c. Approximation C 

Because method C approximates both the streamfunc¬ 
tion and the PV by Galerkin series, the derivation of the 
modal equations is very straightforward compared with 
the calculations in appendix A of Flierl (1978): one sim¬ 
ply substitutes the truncated Galerkin series for the stream- 
function (fT^ and PV (l20l) into the QGPV equation ([fill, and 
then projects onto mode n to obtain 

N N 

“h ^2 P^x^n ~ ^ : (69) 

m=0s=0 

where is defined in (l42l i. and we recall the relation 
between \jf„ and from (l25l l 

qn=^nWn + T,{Pn^^-Pn^^)- ( 60 ) 


In approximation C there are N -|- 3 degrees of freedom: 
the N -h 1 modal amplitudes \j/„ and the two surface buoy¬ 
ancy fields t?^. The approximation C evolution equations 
are completed by advection of the surface buoyancy 

N 

5,t?± + ^p±j(v/„,t?±)=0. (61) 

n=0 

We emphasize that in approximation C the surface buoy¬ 
ancy fields tJ* are not passive: y/n, qn, and t?* are related 
through dbOl l. 

Finally, note that approximation C is recovered from ap¬ 
proximation B if the surface streamfunction is represented 
by a truncated version of the series ( |4^ . 

Quadratic conservation laws 

Approximation C conserves surface buoyancy variance 
as in (Hsi). Total energy is also conserved 

fl\'^Vn\^+l'^n¥ndS = 0. (62) 

“ 11 = 0 -' 

As in approximation B, conservation of potential enstro¬ 
phy is troublesome. With P =0, approximation C has a 
potential enstrophy conservation law 

= (63) 

n=0-' 

But with P ^ 0, approximation C does not conserve the 
analog of the exact potential enstrophy ( fTSl ) (Appendix A). 

5. The Eady problem 

We use classical linear stability problems with nonzero 
surface buoyancy to illustrate how solutions to specific 
problems can be constructed and to assess the relative 
merit and efficiency of approximations A, B, and C. The 
linear analysis does not provide the full picture of con¬ 
vergence of the approximate solutions. Nonetheless, in 
turbulence simulations forced by baroclinic instability, it 
is necessary (but not sufficient) to accurately capture the 
linear stability properties. 

We use nondimensional variables so that the surfaces 
are at = 0 and = —1. The Eady exact base-state 
velocity is given by (l33t with zero PV ^ = 0 and P =0. 

a. Approximation A 

While the surface fields t?^ are dynamically passive in 
approximation A, the Eady problem can still be considered 
because the base-state PV defined via ( l40l i converges to 5- 
distributions on the boundaries (Section^. 

The base-state velocity in Approximation A is given by 
the series ( l34l l and is a good approximation to the exact 
base-state velocity ( l33t . But, according to approximation 
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Fig. 2. Growth rate for the Eady problem as a function of the zonal wavenumber (I = 0) using approximations A, B (exact), C with various number 

of baroclinic modes (N). 


A, there is a nonzero interior base-state PV gradient given 
by the series ( l36l l. As N ^ oo the PV gradient in (l36l l 
converges in a distributional sense to Brethertonian sheets 
at z = 0 and —1. But for numerical implementation of 
approximation A we stop short of N = oo. While the PV 
gradient is much larger at the boundaries, there is always 
interior structure in the PV (Figure [T]:). We show that this 
spurious interior PV gradient has a strong and unpleasant 
effect on the approximate solution of the Eady stability 
problem. 

To solve the Eady linear stability we linearize the 
interior equations (EB about the base-state velocity in 
01 and the PV gradient in (l36l) . We assume = 
qiiex^[\{kx + ly— Ct/'f)], etc, toobtaina (N-l-1) x (N-l-1) 
eigenproblem 


N N 

m =0 j '=0 


Vrr, 


■ Qs ^ 






(64) 


trixform Aq = c^q, where q = ■ ■ • ,?n-i,§n]^ (Ap¬ 

pendix B) and solved with standard methods. 

Eigure |2] shows the growth rate of the Eady instability 
according to approximation A, and compares this with the 
exact Eady growth rate. Approximation A does not do 
well, especially at large wavenumbers. The exact Eady 
growth rate has a high-wavenumber cut-off. At moderate 
values of N, such as 3, 5 and 7 approximation A produces 
unstable “bubbles” of instability at wavenumbers greater 
than the high-wavenumber cut-off The growth rates in 
these bubbles are comparable to the true maximum growth 
rate. As N increases, the unstable bubbles are replaced by 
a long tail of unstable modes with a growth rate that slowly 
increases with k. These spurious high-wavenumber insta¬ 
bilities are due to the rapidly oscillatory interior PV gra¬ 
dient which supports unphysical critical layers: see Eigure 

El 

b. Approximation B, the exact solution 


where Qs are the coefficients of the series ( l36l l and 
jk. The eigenproblem (l64l i can be recast in the ma¬ 


in approximation B, the zero PV in the Eady problem 
implies q„ = ^„= 0. The N-l-1 modal equations (with j3 = 
0) are trivially satisfied; there is no interior contribution 
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Fig. 3. Structure of JC = 8 unstable mode for the Eady problem ob¬ 
tained using approximation A and N = 64. Streamfunction is the black 
curves and PV is the colors. The streamfunction slightly tilts west¬ 
ward as z increases. One can see the unphysical critical layer associated 
with the fast-oscillating base-state PV. The critical level, Zc, is the depth 
where the unstable wave speed matches the velocity of the base-state. 
Only the top quarter of the domain is shown. 


{(j>^ = 0). Thus approximation B solves the Eady problem 
exactly. 

Assuming a = d{z)e^p[i{kx +ly — we obtain 

the solution to the surface streamfunction inversion prob¬ 
lem (|42l)-(|48l) 


' KTsinhif ifsinhK: 


(65) 


where the magnitude of the wavenumber vector is ic = 
\/k'^ + P. We evaluate the surface streamfunction (l65l l at 
the boundaries to find the relationship between the stream- 
function at the surfaces and the boundary fields t?-*; 


c. Approximation C 

Approximation C expands both the streamfunction and 
the PV in standard vertical modes. Thus in the Eady prob¬ 
lem the PV is exactly zero, as it should be: q = qn = 0- 
(This contrasts with approximation A, in which the dif¬ 
ferentiation of the truncated series approximation to the 
streamfunction induces an unphysical oscillatory base- 
state PV gradient.) Thus approximation C does not have 
the spurious critical layers that bedevil A. Moreover, in ap¬ 
proximation C, the N -h 1 modal equations (with j3 = 0) in 
(l59l l are trivially satisfied, and the inversion relationship 
(l60l l provides a simple connection between the stream- 
function and the fields tJ*. The base velocity for the Eady 
problem in approximation C is the series in ( iMl i (the same 
as A). From the exact shear at the boundaries we obtain 
the exact base-state boundary variables 

0± = -y. (70) 

We linearize the boundary equations ( IfiTl i about the base- 
state ( l36l l and dTOl i. to obtain 

, N 

5,tJ± + ^ V4P± = 0. (71) 

k=0 

Assuming tJ* = exp[i{kx + ly — (O^t)], and using the 
inversion relationship (l60l l. we obtain a 2 x 2 eigenprob- 
lem 


’t7+' 

II 

’t9-+' 



p-_ 


where matrix C is defined in appendix C. It is straightfor¬ 
ward to show that converges to the exact eigenspeed. 
i.e., c*' ^ as N ^ oo (Appendix B). Figure|2]shows that 
approximation C successfully captures the structure of the 
Eady growth rate even with modest values of N. 



1 

coth K 

—csch If 


'l>+' 


K 

csch K 

— coth K 


p-_ 


The nondimensional linearized boundary conditions dSTl) 
are 

t?+ - , and -\j/^ = , (67) 

where = afi jk. Using the boundary conditions (l67l i in 

(l66l l we obtain an eigenvalue problem 


1 

K — coth K 

csch K 

'tJ+' 

B 


K 

—csch K 

coth If 

tJ- 

= c 

tJ- 



The eigenvalues of B are given by the celebrated disper¬ 
sion r elation for the Eady problem jPedloskvIl 19871 : IVallisI 
I 2 OO 6 I 1 

= 7±i[(f-^“hf)(f-cothf)] ' . (69) 


d. Remarks on convergence 

The cmdest truncation (i.e. N = 0) is stable for both 
approximations A and C (Figure |2]). With one baroclinic 
mode (N = 1) the growth rates (©,■ = kx Im{c}) are qual¬ 
itatively consistent with the exact solution, and the results 
improve with N = 2. With a moderate number of baro¬ 
clinic modes modes (N > 2) approximations A and C con¬ 
verge rapidly to the exact growth rate at wavenumbers less 
than about 2.2 — see figure |2] But surprisingly the con¬ 
vergence of the growth rate at the most unstable mode 
(if « 1.6) is faster in approximation A ( rxj N-4) than in 
approximation C (~ N^^) — see figure |4] However, the 
convergence in approximation C is uniform: there are no 
spurious high-wavenumber instabilities. 

Figure |4] also shows that the approximation A conver¬ 
gence of the growth rate to zero at if = 8 is slow (~N^’). 
While the growth rate does converge to zero at a fixed 
wavenumber, such as if = 8, we conjecture that there are 
always faster growing modes at larger wavenumbers. 






























































10 


Journal of Physical Oceanography 











Fig. 5. Growth rate for the Green problem with = 1 as a function of the zonal wavenumber (/ = 0) using approximations A, B, C with various 
number of baroclinic modes (N). The black line is a finite-differences solution with 1000 vertical levels. 



Fig. 4. Absolute error as a function of number of baroclinic modes 
(N) for the growth rates of the Eady problem. The solid lines show the 
error at the exact fastest growing mode (fc ss 1.6). The dashed line is the 
approximation A error at K' = 8. 


6. The Green problem 

To further explore the relative merit and efficiency of 
approximations A, B, and C we study the instability prop¬ 
erties of a system with nonzero j3. For simplicity, we con¬ 
sider a problem with Eady’s base-state V'" = — (1 +z)y on 
a g-plane. This i s simil ar to the prob l em or iginally con¬ 
sidered by iGreenI (Il960h and IChamevI lll97 ih . The major 
difference is th at Chamey considered a vertic ally semi¬ 
infinite domain ( ICharney|[l947l : |Pedloskvlll987h while we 
follow Green and consider a finite-depth domain with 
-1 <z<0. 

We obtain the exact system for this “Green problem” by 
linearizing the QG equations ([fill-® about the base-state 
( 1^ with background PV j3y, where 0 is the nondimen- 
sional planetary PV gradient. Assuming yr = 0exp[i(ky-|- 
ly— cot)], we obtain 

{U-c)[$,,-K^^]+0^=O, -l<z<0, (73) 


and 


(t/-c)4-0=O, 


z = -l,0. 


(74) 
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As a reference solution, we solve the eigenproblem 
(I73b - (l74l i using a centered second-order finite-difference 
scheme with 1000 vertical levels: see Figure|5] 

The Green problem supports three classes of unsta¬ 
ble modes, indicated in the lower right panel (N = 128) 
of Figure |5] (1) the “modified Eady modes”, which 

are instabilities that arise from the interaction of Eady- 
like edge waves, only slightly modified by j3; (2) the 
“Green modes”, w hich are very long slowly growing 
modes dVallisIbOO^ : (3) the high-wavenumber “Chamey 
modes” are critical layer instabilities that arise from the in¬ 
teraction of the surface edge wave with the interior Rossby 
wave that is supported by nonzero j3. 

a. Implementation of approximation A 

The base-state for the Green problem is the same as in 
the Eady problem. In approximation A, the j3-term adds 
only a diagonal term to the Eady system (l64l i (see ap¬ 
pendix C). 

b. Implementation of approximation B 

The base-state is the same as in Eady problem. The 
steady streamfunction and buoyancy fields that satisfy dS^ 
and dSTl) exactly are 

E = -(l-|-z)y and 0* = -y. (75) 

Assuming = qn{z)txp[\{kx + ly— CO^t)], the N-f 1 in¬ 
terior equations ( l56l l linearized about ( l75l l become 

N 

^ ^nsqs -f j3 (0„ + a„) =C^qn, (76) 

j=0 

where 

rz'^ 

^ns = S J _PnPs{z+l)dz. (77) 

The boundary conditions dSTl t, linearized about dTSl t, be¬ 
come 

N 

(78) 

,^■=0 

and 

^ (79) 

.^=0 

where d is given by dfiSl) . We use the inversion relation¬ 
ship ( l55l l and the Neumann-to-Dirichlet map ( l66l l to recast 
this eigenproblem into standard form Bq = c^q, where 
q = [B+,qo,qi,...,qN-uq^^,-&-]^ (see appendix C). 

c. Implementation of approximation C 

Again the base-state is the same as in the Eady prob¬ 
lem. But now there are N-f 3 equations: the two boundary 


equations of Eady’s problem ( TTTl i plus N + 1 interior equa¬ 
tions 

N N 

^nms^m ^ ; (80) 

m=05=0 

We use the inversion relationship dfiOl l in dSOl l to recast this 
eigenproblem in the form Cq = c^'q, where q is defined as 
in approximation B (see appendix C). 

d. Remarks on convergence 

The most crude truncation (N = 0) is stable for approx¬ 
imations A and C. In contrast, the N = 0 truncation in ap¬ 
proximation B is qualitatively consistent with the modified 
Eady instabilities: see figure |5] With a moderate number 
of baroclinic modes (N = 2 or 3), approximations A, B 
and C all resolve the modified Eady modes relatively well. 
At the most unstable modified Eady mode (tc « 1.9), ap¬ 
proximation B has typically the smallest error because it 
solves the surface problem exactly. As in the Eady prob¬ 
lem, approximation A converges (~ N^"*) faster than ap¬ 
proximations B and C (~ N^^) at the most unstable mode, 
but B and C converge faster at high wavenumbers (O. 

Approximations A, B, and C all converge very slowly 
to the high-wavenumber Charney modes (Figures |5] and 
These mod es are interior critical-layer instabilities 
( Pedloslc^ Il987h and the critical layer is confined to a 
small region about the steering level (i.e., the depth at 
which the phase speed matches the base velocity — see 
figure |7]i. With finite base-state shear, the critical layer 
is always in the interior. Thus, the problem is not that 
standard vertical modes are inefficient because they do not 
satisfy inhomogeneous boundary conditions; a low resolu¬ 
tion finite-difference solution also presents such “bubbles” 
in high-wavenumber growth rates (not shown). Resolu¬ 
tion of the interior critical layer, not the surface boundary 
condition, is a problem for all metho ds at high wavenum¬ 
bers. The “surface-aware” modes of ISmith and Vannestel 
(l2013h have the same limitation — a large number of ver¬ 
tical modes is required to resolve the critical-layer insta¬ 
bilities (K. S. Smith, pers. comm.). 

For example, with N < 25, at k: = 8, approximations are 
qualitatively inconsistent with the high-resolution finite- 
difference solution. For larger values of N, the growth 
rate convergence for approximations B and C scales ~ 
The growth rate for approximation A converges 
painfully slowly (~N^*). As in the Eady problem, at 
large wavenumbers, the growth rate for approximation A 
is qualitatively different from that of the finite-difference 
solution because of spurious instabilities associated with 
the rapidly oscillatory base-state PV gradient. 

7. Summary and conclusions 

The Galerkin approximations A, B, and C are equiv¬ 
alent if there is are no buoyancy variations at the sur- 
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Fig. 6. Absolute error as a function of number of baroclinic modes 
(N) for the growth rates of the Green problem. The solid line represent 
the error at the exact fastest growing mode (K ^ 1.9). The dashed line 
is the error at /c = 8. 
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Fig. 7. Wave structure of the IC = 8 unstable mode for the Green 
problem with /f = 1 solved using a second-order finite-difference 
scheme with 1000 vertical levels. Streamfunction (black contours) and 
potential vorticity (colors). The streamfunction slightly tilts westward 
as z increases. The potential vorticity is confined to a small region, the 
critical layer. The critical level, Zc, is the depth where the unstable wave 
speed matches the velocity of the base-state. Only the bottom quarter of 
the domain is shown. 


faces. Thus all three approximations are well-suited 


for applications with zero surface buoyancv ( 

Flierllll978l 

IFu and Flierllll980t iHua and Haidvo2ellll986 

). But with 


nonzero surface buoyancy the three approximations are 
fundamentally different. 

App roximation A, originally introduced by iFlierll 
1 197 SI) , represents the streamfunction as a Galerkin se¬ 
ries in standard vertical modes and defines the potential 
vorticity so that the inversion problem is satisfied exactly. 
The most important limitation of A is that the interior PV 
evolves independently of the surface buoyancy (as if t7^ 


were zero). The evolution equations in approximation A 
(HTl i are relatively simple, and the system conserves en¬ 
ergy (|44] | and potential enstrophy (l45l l. But, in approxi¬ 
mation A, nonzero surface buoyancy results in an interior 
PV that distributionally converges to 5-distributions at the 
boundaries as N ^ oo. These smeared-out Bretheitonian 
5-functions provide a very inaccurate representation of 
the true inhomogeneous surface boundary conditi on. Fi¬ 
nite difference schemes have a similar pathology jSmithl 
l2007t Ifulloch and Smithll2009ah . As a result of this arti¬ 
ficial interior PV gradient, solutions with a small number 
of modes are qualitatively misleading and convergence at 
large wavenumbers is very sl ow Slow conver- 

gence was previously noted bv iHua and Haidvogell(ll986h 
— see their figure 2. Furthermore, even if heroic values 
of N achieve convergence at say K = 10, we conjecture 
that there will always be spurious unstable modes at even 
larger wavenumbers. In some simulations these unphysi¬ 
cal high-wavenumber instabilities might be eliminated by 
hyperviscosity or by a scale-selective filter. But one must 
be aware of potential effects on the evolution of the sys¬ 
tem. 

Approximation B, originally introduced by 
iTulloch and Smith! (l2009bh using one baroclinic mode 
and constant buoyancy frequency, takes the opposite 
starting point from approximation A. B represents 
the PV as a Galerkin series in standard modes and 
calculates the streamfunction that satisfies the exact 
inversion problem associated with the approximate PV. 
The linear inversion problem can be split into an interior 
contribution with homogeneous boundary conditions 
and a surface contribution with zero interior source and 
inhomogeneous boundary condi tions (iLanevre and Kleiiil 
l2006t iTulloch and Smith i2009bl) . Thus the exact interior 
streamfunction associated with the approximate PV is 
a Galerkin series, but the surface streamfunction must 
be computed using other methods. Because the surface 
streamfunction projects onto the interior solution the 
energy is not diagonalized. Indeed, approximation B 
conserves neither energy nor potential enstrophy if j3 7 ^ 0. 

Approximation C represents both the PV and the 
streamfunction by Galerkin series. Although the inver¬ 
sion problem is not satisfied exactly, the relation be¬ 
tween modal streamfunction xj/n and PV is obtained by 
Galerkin projection of the exact inversion relationship and 
prominently exhibits the surface buoyancy fields: see dfiOl l. 
Approximation C is the most consistent because it uses the 
same level of approximation for both PV and streamfunc¬ 
tion. The evolution equations ( l59l l are relatively simple, 
and the approximate system conserves total energy (l62l i. 
The most vexing limitation of C is the lack of potential en¬ 
strophy conservation with j3 7 ^ 0. But this does not mean 
that approximation A is better than approximation C: all 
approximations conserve potential enstrophy when j3 = 0, 
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or if T?* = 0 , and none conserve an analog of the exact 
potential enstrophy ddSl) when j3 7 ^ 0 and 7 ^ 0. 

With nonzero interior PV gradients the convergence 
of all approximations is slow for the high-wavenumher 
Charney-type modes. The critical layer associated with 
these modes spans a very small fraction of the total depth 
(Figure |7]i. To accurately resolve these near-singularities 
at the steering level there is no better solution than having 
high vertical resolution in the interior. 

For problems with nonuniform surface buoyancy and 
nonzero interior PV gradient, we recommend approxima¬ 
tion C for obtaining solutions to the three-dimensional QG 
equations using standard vertical modes. 

The codes that produced the numerical re¬ 
sults of this paper, plotting scripts, and sup¬ 
plementary figures are openly available at 
github. com/crocha700/qg_verticaljnodes. 

Acknowledgments. CBR is grateful for a helpful con¬ 
versation with G. R. Flierl. We thank Joe LaCasce 
and Shafer Smith for reviewing this paper. Geoff Vallis 
pointed out that Green (1960) first considered the Eady+j3 
problem. This research was supported by the National Sci¬ 
ence Foundation under award OCE 1357047. 


APPENDIX A 

On the convergence of Galerkin series in standard 
modes 

llacksonl (Il914l) gives conditions for the uniform conver¬ 
gence of series expansions in eigenfunctions of the Sturm- 
Liouville eigenproblem 

^+[p2-A(Z)]P„=0, (Al) 

defined on the interval Z G [0, n] with boundary conditions 

P«(0)-7oP«(0)=0, and K{n) - YnPn{n) = 0 (A2) 

where 70 and 7^1 are real constants of arbitrary sign and 
is the eigenvalue. The equations defining the standard 
modes ([Tli-(|2]i can be brought to this form using the fol¬ 
lowing Liouville transformation 

Z{z) = with 

(A3) 

and 

P„(Z) = 5'(z)‘/'^p„(z), where S{z) = j^^. (A4) 


The boundary condition for the standard modes (O implies 
that the transformed modes satisfy ( IA2b with 


7o = 


As(z-yi^ 

ZdS{z^)/Az'' 


and 'Yk = - 


45(z+)'/2 
ZAS{z+)/Az ' 


(A 6 ) 


If dS '/Az = 0 at a boundary then the appropriate condition 
at that boundary is = 0 . 

A special case of Theorem I from llacksonl (Il914ll states 
that the expansion of a function (j> (Z) as a series in eigen¬ 
functions P„ converges absolutely and uniformly provided 
that both d0/dZ and dA/dZ are continuous and bounded, 
regardless of whether or not (j) satisfies the same bound¬ 
ary conditions as P„. (The remainder of the theorem con¬ 
cerns the rate of convergence under stronger conditions 
on (j) and A.) The streamfunction, potential vorticity, and 
buoyancy frequency profiles are typically assumed to be 
smooth in studies of QG dynamics, which implies that 
both 0 and A will satisfy the above conditions. Uniform 
convergence over Z G [0, ;r] implies uniform convergence 
overz G [z^,z+]. 


APPENDIX B 

Details of the derivation of quadratic conservation 
laws for approximate equations 

a. Approximation A 

To obtain the energy conservation in approximation A 
we multiply the modal equations (|4TI) by — vr„, integrate 
over the horizontal surface, and sum of on n, to obtain 

Xrt [ k[i^Wn? + ^nfn]dS + 

n=0-^ 

N N N . 

J — 0 - (® 1 ) 

tt=0m=05=0 

Notice that 

j ^jtnHVm,Z^,\j/,)AS = j A,vAJ(v>-„,V>m)d5'. (B2) 

Hence the triple sum term in IB II vanishes identically be¬ 
cause is fully symmetric and the Jacobian is skew- 
symmetric. Thus we obtain conservation of energy (l44l) . 
Similarly, to obtain the potential enstrophy conservation 
law in ( 1431 ) we multiply the modal equations (|4TI) by 
AnVA, integrate over the surface, sum on n, and invoke 
the same symmetry arguments used for the energy conser¬ 
vation. 

b. Approximation B 


The eigenvalues are related by p„ = Zk„ and 


A(Z) = Z^ 


1 d^S 
4 Az^ 


1 



Energy nonconservation 


(A5) 


The analog of (fTTIi in approximation B is 
E^=E^+E^+E^^. 


(B3) 
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Eii, = 


The three terms in ( IB 3b are 

E /1 l\n 

n=0 


dV 




d5, 


E„ = 


i/i l^<Jl" + (^) M 


O) 


dV, 


(B4) 


(B5) 


and 


E^o = };J V</.«.Va+(|) d,<j4d,a 

N 


dV, 


,d5. 


IN /* 

- L L i /p«P..(0« + O'«)J(c^,?.OdV 

n=0s=0 

N N N „ 

+ L LL ^mns I J(0m; dS'. 


df 


a = 


cosh (z+ 1 ) 


COSA- 


coshz 


■ sinA. 


7^0. 


(B 6 ) 


The cross-term £^<7 is not zero because the surface stream- 
function a projects on the standard vertical modes i.e., 0 ,, 
is nonzero. To obtain an equation for Eb we form evolu¬ 
tion equations for the three components in ( IB3b and add 
them. The final result is 
dE^ N 


POTENTIAL ENSTROPHY NONCONSERVATION 

The analog of the exact potential enstrophy (fTSl l is not 
conserved in approximation B. We attempt to form a po¬ 
tential enstrophy conservation by multiplying the interior 
equations by A„ 0 „ and integrating over the surface, 
and summing on n, 

d ^ r . ^ /■ 

TtL / ji^n^nfdS - P Y, / A„d„dj,^„dS = 0. 

°'n=0"' 11=0-' 

(Bll) 

The potential enstrophy given by the sum on the left-hand- 
side of (IB 11b is conserved in the special case P = 0. For 
j3 7 ^ 0 we first form an equation for A„a„, and then cross- 
multiply with the modal equations dSfil l. integrate over the 
surface, and sum on n. Eliminating the j3-term in (IB 11b 
yields 


df 


Y + {A„6„)qndS = 


n=0 


N N N 


(B7) 


=z0n=0s—0 

The simplest model with barotropic interior dynamics 
(N = 0) conserves energy. With richer interior stmcture, 
however, the right-hand-side of (IB7b is generally nonzero. 
Consider the “two surfaces and tw o modes” (TMTS) 
model of iTulloch and SmithI ( l2009bl) . corresponding to 
N = 1 with constant buoyancy frequency. Using nondi- 
mensional variables the energy equation (IB7b becomes 

d£f 


-f 


IN IN IN p 

.n=0n=0s=0 

F7=0 


«=0m=0 


-p„J(a -fp„,0„,,tJ )]d5. 

(B12) 


y 1^01 J (ai, Ai0i) -J (cti , 0 - 2 )] dS. 

(B 8 ) 

We now construct an example in which we can analytically 
show that the right-hand-side of (IB 8 b is nonzero. This ex¬ 
ample should be interpret as an initial condition for which 
energy is guaranteed to grow or decay. For simplicity we 
consider A 0 i = A. (j>i, where A is a constant, so that the first 
term on the right-hand-side of (IB 8 b is identically zero. As 
for the surface streamfunction, we choose 


The right-hand-side of (IB 12b is nonzero even in the sim¬ 
plest model (N = 0). 

c. Approximation C 

To obtain the conservation of energy in approximation 
C we multiply the modal equations ( l59l ) by — V/,,, integrate 
over the horizontal surface, and sum on n, to obtain 


A "" 

df 


(B9) 


sinh 1 sinh 1 

We use 01 = sin A cosy so that all fields are periodic with 
the same period. Integrating over one period in both direc¬ 
tions gives 

d£f A 

~ V2{\+5n'^ + An^) 

The total energy Ef grows or decays depending on the 
sign of A. Thus the analog of the exact energy (fTTb is not 
generally conserved in approximation B. 


(BIO) 


Y [ [i^Wnf + ^^f^„]dS 

n=0'^ 

-Hi [ fndtip+^+-p;;^^)ds 

«=o 

N N N . 

«=0m=0j=0 

The triple sum term vanishes by the same symmetry ar¬ 
guments used above in approximation A. The term on the 
second line of (IB 13b is also zero; multiply the boundary 
conditions dMT l by p* ij/n and integrate over the horizontal 
surface. Thus we obtain the energy conservation law in 
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Potential enstrophy nonconservation 


c. Approximation B 


As in approximation B, the analog of the exact potential 
enstrophy (fTSl l is not conserved in approximation C. The 
potential enstrophy equation with j3 7 ^ 0 is formed anal¬ 
ogously to the approach used above in approximation B. 
The final result is 

d N ^ -2 

-f+q„Ana„dS = 

N N . 

+ L Li ^nPnPV{f'n,'&+) 

m = Qn =0 
N N N . 

“ / ^mns^n^n-^ i^s) ■ (B14) 

m=0n=0 5=0'^ 

The right-hand-side of (IB 14b is zero for the simplest 
model (N = 0), but it is generally nonzero. 


APPENDIX C 

Details of the stability problems 

a. The interaction tensor 

Because the standard vertical modes with constant strat¬ 
ification are simple sinusoids (l32l i. the interaction coeffi¬ 
cients (I 42 I 1 can be computed analytically. First we recall 
that Sjjj. is fully symmetric. Permuting the indices so that 
i > j > k 'we obtain 


^ijk = < 


1 ; 

2 ' 

0 : 


i = j,k = 0; 
i = j + k- 
otherwise. 


(Cl) 


The second line in tClb c orrects a factor of i missed by 
iHua and Haidvog^ ( 19861) . 


The Green eigenvalue problem in d?^ through ( l79l ) 
can be recast in the standard form Bq = c^q, where 
q = ,qo,qi ,... The first and last rows 

of the system stem from the boundary conditions (l78l l-(l79]l 


, cothifX ^ cschic - 

1 - —— t?+ - Y^ptaS - —— B- 

/C / /c 

/ 5=0 




(C4) 


and 


cschif - , ^ _ cothK: 

-t?+-^p, - B . (C5) 

^ .5=0 ^ 

The {n + l)’th row originates from the «’th interior 
equation (|76| 


^ Pn Yj Ynsqs “f {l^^n ^ 1) 4“ ^ P« ^ — C qn i 

5=0 

(C 6 ) 


where the symmetric matrix y,ns is 


i . 

■ 2 • 




P,P,zdz= TJW'- ' = 0-4isodd; 


i + j is odd. 


(C7) 


d. Approximation C 
The Eady problem 
The 2x2 eigenproblem is 


+En — 

'tJ+' 


't?+' 



— L 



(C 8 ) 


b. Approximation A 

Using the symmetry in and the inversion relation 
(HSll, we rewrite row n + 1 of the linear Green system 

N N 

^nms {Um -|- B\Qm 0!^^ C^s H“ ^GCn^n — ^Qni (C2) 

5=0m=0 

where the inverse of the «’th mode Helmholtz operator in 
Fourier space is 


a„ =-(K-^-l-(«;r)^) ‘. (C3) 

The Eady problem is the special case (5=0. We use a 
standard eigenvalue-eigenvector algorithm to obtain the 
approximate eigenspeed . 


where 

. N . N 

EN = ao + 2 y^a„, and = ao + 2 ^ (-l)'’a„. 

H=1 n =\ 


(C9) 


For finite N, the approximate eigenspeed is 


r/G+ I jjG 

^ + On _[_ 


r/G+ tjG \ 2 




1/2 


. (CIO) 


The sums (IC9b become exact in the limit N ^ oo 


E« = - 


cothic 


K 


and n. = -£!£!!£. (cii) 
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The base velocity also converges to the exact result. No¬ 
tice that 


11 f 1 

+ + ~ 8 

and therefore 

U2^ = l, and C/°^=0. 


(C12) 


(C13) 


Thus 

C^-B as N^oo, (C14) 

and the eigenvalues of the Eady problem using approxi¬ 
mation C become exact i.e., c*' ^ as N ^ oo. 


The Green problem 
The (N -P 3) X (N -P 3) eigenproblem is 

Cq=c‘'q, (C15) 

where q is defined as above in approximation B. The first 
and last rows of ( IClSb stem from the boundary conditions 

dlD 


(u2^+L^'^ -'^a„p2q„-Q.N^ =c^i^+, 

n=0 

(C16) 

and 

N /-\- 

a„p^^„-P - Ln j . 

n=0 

(C17) 

Row «-P 1 originates from the «’th modal equation dSOl i: 

N N 

5=0m=0 

—j3a„p^ t?^ = c‘'^„. (C18) 
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